$Title  CF_CO2_NH.GMS

* ###################################

* Per capita income, consumption patterns, and $CO_2$ emissions

* Journal of the Association of Environmental and Resource Economists

* Justin Caron and Thibault Fally

* April 2021

*###################################


* COMPUTES FITTED PRODUCTION SHARES (consistent with fitted consumption shares) FOR USE IN SIMULATION MODEL


$set SLASH \
$if %system.filesys% == UNIX $set SLASH /


* declare set of demand estimates to use : theta4, notc, tc

$if not set demandest $set demandest tc
* need to choose H or NH here


* CRIE: NH
*$if not set spec $set spec NH

* Or CLM:

* options: NH_CLM
* options: NH_CLMquad
* options: NH_CLMshift

$if not set spec $set spec NH
$if not set ds $set ds gtap8gas

* TWO COUNTERFACTUALS: CHOOSE ONE:
parameters income_shock, tcost_shock, income_shock_toresourcefactor;
income_shock = 1;
tcost_shock = 1;
income_shock = 1.01;

* set this to 1 if productivity of resource factor is to be shocked as well -- 0 if off
income_shock_toresourcefactor = 1;
$if "%prodspec%"=="noresprod"  income_shock_toresourcefactor = 0;

* uses data uploaded in data_CF_Nov2015.gdx
* variables:
* sigma, theta, beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm;

Sets
r country  /
AUS,NZL,CHN,HKG,JPN,KOR,MNG,TWN,KHM,IDN,LAO,MYS,PHL,SGP,THA,VNM,BGD,IND,NPL,PAK,LKA,CAN,
USA,MEX,ARG,BOL,BRA,CHL,COL,ECU,PRY,PER,URY,VEN,CRI,GTM,HND,NIC,PAN,SLV,AUT,BEL,CYP,CZE,
DNK,EST,FIN,FRA,DEU,GRC,HUN,IRL,ITA,LVA,LTU,LUX,MLT,NLD,POL,PRT,SVK,SVN,ESP,SWE,GBR,CHE,
NOR,ALB,BGR,BLR,HRV,ROU,RUS,UKR,KAZ,KGZ,ARM,AZE,GEO,BHR,IRN,KWT,QAT,SAU,TUR,ARE,EGY,MAR,
TUN,CMR,CIV,GHA,NGA,SEN,ETH,KEN,MDG,MWI,MUS,MOZ,TZA,UGA,ZMB,ZWE,BWA,NAM,ZAF,ISR,OMN
/;


alias (r,s);
alias (r,rloop);

Sets
i industry  /
isr,obs,ros,osg,pdr,wht,gro,v_f,osd,c_b,pfb,ocr,ctl,oap,rmk,wol,frs,fsh,coa,oil,gas,omn,cmt,omt,vol,mil,pcr,
sgr,ofd,b_t,tex,wap,lea,lum,ppp,p_c,crp,nmm,i_s,nfm,fmp,mvh,otn,ele,ome,omf,ely,wtr,cns,trd,otp,wtp,atp,cmn,ofi
/;

alias (i,j);

Sets
f factor /Land,UnskLab,SkLab,Capi	tal,NatlRes/;

* Loading data:
parameter fittedexp, coeffs, scale, sigma, theta, beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop, incelast;

* Load demand estimates from DEMAND.GMS
* CRIE
$if "%spec%"=="NH" $gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_rall.gdx


* If loading CLM
$if "%spec%"=="NH_CLMshift"  $gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_CLM_flex_shifter.gdx
$if "%spec%"=="NH_CLMquad"  $gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_CLM_flex_quadratic.gdx
$if "%spec%"=="NH_CLM"  $gdxin estimates%SLASH%estimates_%ds%_logweighted_%demandest%_CLM.gdx

$load  coeffs scale= sigma_scale fittedexp

* And everything else from DATAPREPARATION.GMS

$gdxin estimates%SLASH%co2_data_%demandest%.gdx
$load  beta, intersh, trade, endow, vom, vdfm, vifm, fd, vfm, pcexp, pop, incelast

theta(i) = 4;

* -------
* Theta 4:  H, sigma = scale, NH, sigma = sigma * scale
*  here, the average sigma level is recoverd from Mu

* TC: H, we don't know the average level of sigma
* back it out assuming the average level of theta = 4

* IF USING BACKED OUT THETA:
*$if "%demandest%"=="tc" theta(i)$coeffs("backed out theta","coeff",i) = coeffs("backed out theta","coeff",i);

* IF USING THETA = 4 IN THE SIMULATIONS
$if "%demandest%"=="tc" theta(i) = 4;


*normalize theta's such that the average equals = 4 and min = 1:
theta(i)$(theta(i) = 0)  = 0;
theta(i)$(theta(i) < 0)  = 0;
theta(i) = 4 * theta(i) * sum(j,1) / sum(j,theta(j));
theta(i)$(theta(i) = 1)  = 1;
theta(i)$(theta(i) < 1)  = 1;
theta(i) = 4 * theta(i) * sum(j,1) / sum(j,theta(j));


theta(i)$(theta(i) > 12)  = 12;
* no change here


* TC - RENORMALIZE SO THAT WE GET CORRECT AVERAGE PRICE ELAST 

sigma(i) =  coeffs("sigma","coeff",i);


parameter  sigma_rescale;
sigma_rescale = (sum(j$sigma(j), 1) + 4 * sum(j, coeffs("mu","coeff",j) )) / sum(j$sigma(j), sigma(j));
display sigma_rescale;

* TC
* rescale sigmas such that average price elasticity (1-sigma) = average mu * theta
* leads to high sigmas
$if "%demandest%"=="tc" $if "%spec%"=="NH"  sigma(i) =  coeffs("sigma","coeff",i);
$if "%demandest%"=="tc" $if "%spec%"=="NH" sigma(i) = sigma_rescale * coeffs("sigma","coeff",i)


* here, in H, use the same average sigma (could also recompute using estimated mu's with H prefs..)
$if "%demandest%"=="tc" $if "%spec%"=="H"  sigma(i) = sigma_rescale;


$if "%demandest%"=="theta4" $if "%spec%"=="NH"  sigma(i) = scale("non-homoth") * coeffs("sigma","coeff",i);
$if "%demandest%"=="theta8" $if "%spec%"=="NH"  sigma(i) = scale("non-homoth") * coeffs("sigma","coeff",i);


*T4: sigma = 2.8
$if "%demandest%"=="theta4" $if "%spec%"=="H"  sigma(i) = scale("homoth") ;
$if "%demandest%"=="theta8" $if "%spec%"=="H"  sigma(i) = scale("homoth") ;


set not_fd(i)  sectors with no sigma estimate (where dropped from estimation);
not_fd(i)$(not sigma(i))  = yes;


* if no sigma, assign average sigma
sigma(i)$not_fd(i)  = sum(i.local$sigma(i), sigma(i)) / sum(i.local$sigma(i), 1) ;

* set theta = 4 for sectors for which we have no estimates
theta(i)$(not theta(i)) = 4;

* or simply impose value 4 to all theta's:
display theta, sigma;


* ----
* CHOICE OF PRODUCTION SHARES 

* compute fitted production shares from fitted demand shares

* 1- load fitted demand vector
* Loading data:
set fit /non-homoth, homoth, true/;


** REDEFINING ABSORPTION AND PRODUCTION FOR CONSISTENCY:
parameters production, absorption;

production(i,r) = sum(s, trade("data",i,r,s));

absorption(j,r) = fd(j,r) + sum(i,intersh(j,i,r) * production(i,r));

parameter
dem_sh(i,r)  share of sector i in total final demand
trade_sh(i,r,s)  share of country r in country s's absorption (Import shares)
prod_sh(i,r,s)  share of country s in country r's sales
fact_sh(f,i,r)   share of sector i in demand for factor f
endow_sh(f,r) share of factor f in total endowment (income)
gosh_m(j,i,r) gosh matrix coefficient (how much j is used in industry i)
gosh_f(j,r) gosh matrix coefficient (how much j is used by final consumers);


* j being used in industry i:
gosh_m(j,i,r)$absorption(j,r) = intersh(j,i,r) * production(i,r) / absorption(j,r);

gosh_f(j,r) = 1 - sum(i,gosh_m(j,i,r));
* normalized to unity when absorption is zero

* check that there is no negative value:
display gosh_f;
* Some negative values because of rounding errors, not a problem 


dem_sh(i,r) =  fd(i,r) / sum(i.local, fd(i,r));

* defined as: from r to s:
trade_sh(i,r,s)$sum(r.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(r.local , trade("data",i,r,s));

parameter comp_ab1, comp_ab2;
comp_ab1(i,s)$absorption(i,s) = round(sum(r,trade_sh(i,r,s)) - 1,14);
comp_ab2(i,s)$(round(sum(r,trade_sh(i,r,s)) - 1,14)) = absorption(i,s);
display comp_ab1, comp_ab2;


trade_sh(i,s,s)$(sum(r.local,trade("data",i,r,s))=0)  = 1;

* defined as: from r to s:
prod_sh(i,r,s)$sum(s.local , trade("data",i,r,s))  = trade("data",i,r,s) / sum(s.local , trade("data",i,r,s));


* ADJUSTMENT WHEN SUM IS ZERO:
prod_sh(i,r,r)$(sum(s.local , trade("data",i,r,s))=0)  = 1;

fact_sh(f,i,r)$sum(i.local,  beta(f,i,r) * vom(i,r)) =  beta(f,i, r) * vom(i,r) / sum(i.local,  beta(f,i,r) * vom(i,r))  ;

endow_sh(f,r) = sum(i,vfm(f,i,r)) / sum((f.local,i),vfm(f,i,r)) ;

* CHECKS TO ENSURE THAT CHG=1 IS A SOLUTION:

parameters CHECK_budget, CHECK_absorption, CHECK_production, CHECK_supplier, CHECK_index, CHECK_fmc, CHECK_income;
CHECK_budget(r)= round(sum(i,dem_sh(i,r)) - 1,14);
CHECK_absorption(j,r) = round(gosh_f(j,r) + sum(i$gosh_m(j,i,r),gosh_m(j,i,r)) - 1,14);
CHECK_production(i,r) = round(sum(s,prod_sh(i,r,s)) - 1,14);
CHECK_supplier(i,r) = round(prod(f$beta(f,i,r),1**beta(f,i,r))**theta(i) * prod(j,1**intersh(j,i,r))**theta(i) - 1,14);
CHECK_index(i,s) = round(sum(r,trade_sh(i,r,s)) - 1,14);
CHECK_fmc(f,r) = round(sum(i,fact_sh(f,i,r)) - 1,14);
CHECK_income(r) = round(sum(f,endow_sh(f,r)) - 1,14);

display CHECK_budget, CHECK_absorption, CHECK_production, CHECK_supplier, CHECK_index, CHECK_fmc, CHECK_income;

* Now, recreating absorption and production given fitted (final) demand and observed import shares:
* quick model to solve for production, based on fitted consumption and observed trade shares and IO links:
positive variables reconstruct_absorption(i,r) fitted absorption
                   reconstruct_production(i,r) fitted production;

parameter final_demand, trade_sh_adj;
final_demand(i,r) = fd(i,r);
trade_sh_adj(i,r,s)$absorption(i,s) = trade("data",i,r,s)/absorption(i,s);


equations eq_IOreconstruct, eq_prodreconstruct;

eq_prodreconstruct(i,r)..   reconstruct_production(i,r) =e= sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption(i,s) );
eq_IOreconstruct(i,r)..     reconstruct_absorption(i,r) =e= final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production(j,r));


model Reconstruct / eq_prodreconstruct.reconstruct_production, eq_IOreconstruct.reconstruct_absorption /;
Reconstruct.iterlim = 1000000;

Reconstruct.reslim = 1000000;



reconstruct_production.L(i,r) = sum(s.local , trade("data",i,r,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));

solve Reconstruct using MCP ;

parameters CHECK_production_bis;

CHECK_production_bis(i,r)$reconstruct_production.L(i,r) = round(1 - sum(s.local , trade("data",i,r,s)) / reconstruct_production.L(i,r), 14);
display CHECK_production_bis;


parameters fittedprod, fitted_prod_sh;

fittedprod(i,r,"true") = reconstruct_production.L(i,r);
fitted_prod_sh(i,r,s,"true")$reconstruct_production.L(i,r) = trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s) / (sum(s.local, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s)));
fitted_prod_sh(i,r,r,"true")$(reconstruct_production.L(i,r)=0) = 1;

* not needed
parameter fitted_production_absorption;
fitted_production_absorption("absorption","true",i,s) = reconstruct_absorption.L(i,s);
fitted_production_absorption("production","true",i,r) = reconstruct_production.L(i,r);


* NOW: COMPUTING  EXPORT AND PRODUCTION PATTERNS with NH prefs:
* note: problems when fitted final demand is nonzero but true final demand is zero:

final_demand(i,r)$fd(i,r) = fittedexp("non-homoth",i,r);


reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));

solve Reconstruct using MCP ;

fittedprod(i,r,"non-homoth") = reconstruct_production.L(i,r);
fitted_prod_sh(i,r,s,"non-homoth")$reconstruct_production.L(i,r) = trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s) / (sum(s.local, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s)));
fitted_prod_sh(i,r,r,"non-homoth")$(reconstruct_production.L(i,r)=0) = 1;

fitted_production_absorption("absorption","non-homoth",i,s) = reconstruct_absorption.L(i,s);
fitted_production_absorption("production","non-homoth",i,r) = reconstruct_production.L(i,r);


parameter distortion_prod;
distortion_prod(i,r,"non-homoth")$fittedprod(i,r,"true") = fittedprod(i,r,"non-homoth") / fittedprod(i,r,"true");

DISPLAY distortion_prod;


* NOW: COMPUTING  EXPORT AND PRODUCTION PATTERNS with H prefs:
* skip for CLM:

$if "%spec%"=="NH_CLMshift" $goto skip
$if "%spec%"=="NH_CLMquad"  $goto skip
$if "%spec%"=="NH_CLM"  $goto skip

final_demand(i,r)$fd(i,r) = fittedexp("homoth",i,r);


reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));
reconstruct_production.L(i,r) = sum(s, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s));
reconstruct_absorption.L(i,r) = final_demand(i,r) + sum(j, intersh(i,j,r)*reconstruct_production.L(j,r));

solve Reconstruct using MCP ;

fittedprod(i,r,"homoth") = reconstruct_production.L(i,r);
fitted_prod_sh(i,r,s,"homoth")$reconstruct_production.L(i,r) = trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s) / (sum(s.local, trade_sh_adj(i,r,s)*reconstruct_absorption.L(i,s)));
fitted_prod_sh(i,r,r,"homoth")$(reconstruct_production.L(i,r)=0) = 1;

fitted_production_absorption("absorption","homoth",i,s) = reconstruct_absorption.L(i,s);
fitted_production_absorption("production","homoth",i,r) = reconstruct_production.L(i,r);


parameter distortion_prod;
distortion_prod(i,r,"homoth")$fittedprod(i,r,"true") = fittedprod(i,r,"homoth") / fittedprod(i,r,"true");

DISPLAY distortion_prod;

$label skip


execute_unload 'estimates%SLASH%fitted_prod_shares_%spec%_%demandest%.gdx',
fitted_prod_sh, fittedexp, distortion_prod, fitted_production_absorption, final_demand	;




